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ABSTRACT 



The formation of the first stars out of metal-free gas appears to result in stars 
at least an order of magnitude more massive than in the present-day case. We 
here consider what controls the transition from a primordial to a modern initial 
mass function. It has been proposed that this occurs when effective metal line 
cooling occurs at a metallicity threshold of Z/Zq > 10~^'^. We study the influence 
of low levels of metal enrichment on the cooling and collapse of initially ionized 
gas in small protogalactic halos using three-dimensional, smoothed particle hy- 
drodynamics simulations with particle splitting. Our initial conditions represent 
protogalaxies forming within a previously ionized H ll region that has not yet had 
time to cool and recombine. These differ considerably from those used in sim- 
ulations predicting a metallicity threshold, where the gas was initially cold and 
only partially ionized. In the centrally condensed potential that we study here, 
a wide variety of initial conditions for the gas yield a monolithic central collapse. 
Our models show no fragmentation during collapse to number densities as high 
as 10^ cm~^, for metallicities reaching as high as 10~^ Zq, far above the threshold 
suggested by previous work. Rotation allows for the formation of gravitationally 
stable gas disks over large fractions of the local Hubble time. Turbulence slows 
the growth of the central density slightly, but both spherically symmetric and 
turbulent initial conditions collapse and form a single sink particle. We therefore 
argue that fragmentation at moderate density depends on the initial conditions 
for star formation more than on the metal abundances present. The actual initial 
conditions to be considered still need to be determined in detail by observation 
and modeling of galaxy formation. Metal abundance may still drive fragmen- 
tation at very high densities due to dust cooling, perhaps giving an alternative 
metallicity threshold. 

Subject headings: stars: formation - stars: mass function - early universe - 
hydrodynamics - equation of state - methods: numerical 



Introduction 



Observations of old stellar populat ions reveal no primordial, metal-free stars fiBeers fc Christlieb 



20051 : iFrebel. Johnson fc Bromnj|2007l ) . although the lowest-metallicity stars have abundance 
distributions suggesting t hat they may have formed from gas polluted by the ejecta of a 
single supernova (see e.g. iTominaga. Umeda fc Nomotd 120071 ). Models of primordial star 
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formation appear to have converged on the resuh that prim ordial stars formed w i th an 



Abel. Brvan. fc Norman 



20061: 



2002 



dern stars ( 


Omukai & Palla 


2001. 


2003: 


O'Shea & Norman 


2006: 


Yoshida et al. 



McKee fc Tarul2008l ). Rather than predominantly forming stars with masses similar to 



or less than the Sun's, primordial stars seem to have been massive, short-lived, objects. 

This raises t he question o f what controls the transition from the primordial IMF to 
the modern one. iBromm et al.l (120011 ) performed simulations of the collapse of cold gas in 
a top-hat potential that included the metallicity-dependent effects of atomic fine structure 
cooling. In the absence of molecular cooling, they found that fragmentation suggestive 
of a modern IMF only set in at metallicities above a threshold value of Zth — lQ~^Zp) . 



Howev er, they noted that the neglect of molecular cooling could be significant. lOmukai et al. 



(120051 ) argued, based on the results of their detailed one- zone models, that mo lecular cooling 



would i ndeed dominate the cooling over many orders of magnitude in density. IJappsen et al. 



(j2007bl . hereafter Paper IV) presented the results of three-dimensional collapse simulations 
that included molecular cooling modeled using th e simpli f ied ch emical network described in 



over 



J2OO8I. le reafte r Paper III). They 



Bromm et al.l (120011 ) occurs in models 



Glover &: Jappsen] ( 120071 . hereafter Paper I) and [G 
found that fragmentation similar to that seen by 
starting from the same initial conditions, but with metallicities below the threshold, and 
indeed even with zero metallicity. 



These results suggest that the initial conditions adopted by iBromm et al.l (I2OOII ) may 
have determined the result much more than it might have been appreciated at the time. 
These initial conditions can now be seen to represent not necessarily the predominant mode 
of early low-metallicity star formation: the gas begins at a redshift oi z = 100 with its 
metallicity already set, but with its temperature the same as that of unenriched gas at 
that redshift, and far below the virial temperature of the halo into which it is eventually 
incorporated. Furthermore, the dark matter halo itself has no central condensation, but 
merely local perturbations on a fiat, top-hat, potential. 

It is therefore important to consider whether a metallicity threshold appears iii sim- 
ulations with different initial conditions than those adopted by iBromm et al.l (I2OOII ) and 
Paper IV. In this paper, we examine a different, still idealized, set of initial conditions, and 
find no fragmentation even at metallicities above the claimed threshold. We therefore argue 
that fragmentation at moderate density depends on the initial conditions for star formation 
more than on the metal abundances present. The actual initial conditions to be considered 
still need to be determined in detail by observation and modeling of galaxy for mation. Metal 



abundance may still dri v e frag mentation at very high densities as suggested by lOmukai et al. 



( I2OO5I ) and IClark et al.l (120081 ). perhaps giving an alternative metallicity threshold. 
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In ^we describe the numerical methods and initial conditions used for our models, while 
in ||3]we describe our results. In [|l]we compare our results to earlier work to demonstrate 
the importance of the environment and initial conditions to the problem. 



Numerical Methods 



2.1. Algorithms 



To assess the relative influence of initial conditions and metallicity on the cooling and 
collapse of gas in small protogalactic halos, we used numerical simulations. During collapse, 
gas increases in density by several orders of magnitude, and so it is best simulated by 
a numerical method with a high dynamical range. We therefore chose smoothed particle 
hydrodynamics (SPH). Overviews of the r netho d , its n ume rical implementati on, and some 
of its applications are given in reviews b y iBeng (119901) and iMonaghanI (119921 ). We use the 
parallel SPH code GADGET version 1.1 (jSpringel et al.ll200ll ) for our simulations. 



SPH is a Lagrangian method for simulating astrophysical flows, in which the fluid is 
represented by an ensemble of particles, with flow quantities at a particular point obtained 
by averaging over an appropriate subset of neighboring SPH particles. The mass resolution 
of a simulation is approximately 

Mres = 2Aneigh'^p (1) 

where N^eigh is the number of particles within a given SPH smoothing kernel and rrip the 
mass of a single gas particle. To avoid numerical fragmentation the Jeans mass Mj has to be 
resolved (iBate &: Burkertlll997l ): Mres < Mj. In order to achieve a higher mass resolution, 
we refine t he mass of the gas particles in regions approaching the Jeans criterion using the 
method of iKitsionas fc WhitworthI (120021 ). For details on the implementation of the method 
see Appendix lAl 

Particle timesteps are constrained by the Courant-Friedrich-Lewy condition that signals 
cross no more than a fraction of a resolution element per timestep, so they grow increasingly 
short as the resolution improves during collapse of high density regions. Replacing dense 
cores with artificial sink particles therefore leads to a considerable increase in computational 
speed, allowing the dynamical evolution of the lower-density gas to be followed over multiple 
free-fall times. In the runs pre sented here we implement sink particles according to the 
prescription of iBate et al.l (119951 ). Sink particles are introduced in regions where the density 
rises above 1.25 x 10^ cm~^, and accrete gas within a radius of 0.1 pc. On every timestep, any 
gas within this accretion radius that is gravitationally bound to the sink particle is accreted 
by it. The design and implementation of our sink particle algorithm is discussed in more 
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detail in lJappsen et al.l (120051 ) . 



2.2. Chemistry and Cooling 

2.2.1. Chemical Model 

The chemical model used in our simulations is the same as that used in paper IV, and 
is discussed in more detail in papers I & III. It is designed to accurately follow the major 
atomic and molecular coolants of the gas. Provided that carbon and oxygen are amongst 
the most abundant metals, the major coolants will be largely the same as i n local atomic 



and molecular gas, namely H2, HD, C, C"*", O, Si, Si"*", CO, OH and H2O (lOmukai et al. 



20051 ) . We therefore follow the abundances of these ten species, together with an additional 
29 species that play important roles in determining the abundances of one or more of these 
coolants. A full list of the chemical species included is given in Table [H The chemical 
network presented in papers I and III contains 189 collisional gas-phase reactions between 
these 39 species, as well as grain surface reactions, reactions involving the photoionization or 
photodissociation of chemical species by ultraviolet radiation, and reactions involving cosmic 
rays. For simplicity, in the simulations presented in this paper we do not include the effects 
of dust, UV radiation or cosmic rays, and so use a simplified version of the model that 
contains only the collisional reactions. 

To implement this chemical network within GADGET, we make use of operator splitting. 
During a given SPH particle time-step, we first compute the new SPH densities in the 
standard fashion, and then update the chemical abundances by solving a coupled set of rate 
equations of the form 

^ = Q- D,n., (2) 
at 

where rii is the number density of species i, and Ci and Di are chemical creation and de- 
struction terms that generally depend on the temperature T and the chemical abundances 
of the other reactants in the system. In our current implementation, we solve rate equations 
for the abundances of only 18 of our 39 species. The abundances of a further 14 species 
that have rapid creation and destruction timescales are determined under the assumption 
that chemical equilibrium applies, while the abundances of the final seven species are deter- 
mined through the use of conservation laws for charge and elemental abundance. For each 
chemical species in our model, the method of solution used is summarized in Table [H To 
ensure numerica l stability, we solve our set of chem ical rate equations implicitly, using the 



DVODE solver (IBrown. Byrne, fc Hindmarshlll989l ). together with an implicit equation for 



the specific internal energy of the particle (see §2.2.21 below). 
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2.2.2. Cooling Function 

The thermal evolution of the gas in our simulations is modelled using a cooling function 
that includes the effects of atomic fine structure cooling from C, C"*", O, Si and Si^, rotational 
and vibrational cooling from H2, HD, CO and H2O, Lyman-a cooling, Compton cooling, and 
H+ recombination cooling, as well as a number of other processes of lesser importance. A 
full list of the processes included is given in Table [21 further details can be found in papers I 
and III. To allow for the effects of cosmic microwave background (CMB) heating of the gas, 
we adopt in our simulations a modified cooling rate of the form 

A = A(r)-A(TcMB), (3) 

where A(T) is the cooling rate of the gas per unit volume for gas temperature T and A(Tcmb) 
is the cooling rate when T = Tcmb- This modification has a negligible effect when T ^ Tcmb, 
but prevents the gas from cooling beneath the floor set by the CMB temperature, except by 
adiabatic expansion. 

To treat radi ative cooling within the GADGET framework, we use the same isochoric 



approximation as ISpringel et al.l (l200ll ). During a given particle time-step, we first compute 
Mad, the rate of change of the internal energy due to adiabatic gas physics. We then solve an 
implicit equation for the new internal energy: 



A [p",m"+^] At 



u-+^ = + u^M - ' ^ , (4) 

pn 

where m" and u^^^ are the internal energy per unit mass at time t" and t""^^ respectively, 
p" is the gas density at time t". This implicit equation is solved simultaneously with the 
chemical rate equations using the DVODE solver. 

We ensure that the internal energy of the SPH particle does not change by a large 
amount during a single timestep by constraining the timestep: 

At < r/c X mm — , — — , (5) 

V^^ad A / 

where rjc = 0.01. This limitation is necessary to ensure that the estimates constructed by 
GADGET of the internal energy and thermal pressure at intermediate points during the 
timestep are accurate, and hence that the pressure forces acting on particles with shorter 
timesteps are also computed accurately. 
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2.3. Initial Conditions 



2.3.1. Initial Temperature and Density Distribution 



Our initial conditions are based on those used in IJappsen et al.l (l2007al hereafter Paper 
11), although collapse is followed to much higher density in the current work. We study 
protogalaxies forming from gas with varying degree of metallicity. We assume the gas is 
initially fully ionized, with no molecules present, and with te mperature Tg = 10'^ K. These 



initial conditions are intended to represent a fossil H ll region (jOh fc Haimanll2003l ) polluted 
by supernovae from the ionizing object. Because of the higher mixing rates expected in gas 
with higher sound speed, such warm gas will have metals mixed into it far earlier than cold 
primordial gas in the same region. Fossil H ll regions subject to cooling and collapse will be 
common as the characteristic lifetimes of ionizing sources are significantly shorter than the 
Hubble time even at high redshifts. The initial temperature of gas in such a fossil H ll region 
should, strictly speaking, depend on various factors, such as the spectrum and strength of 
the ionizing source and the density of the gas. However, since the gas cools rapidly through 
Ly-a cooling at the start of our simulations and since the temperature at which Ly-a cooling 
becomes ineffective does not depend on the initial temperature of the gas, we do not expect 
our results to be sensitive to small changes in the initial value of Tg. Because of the cooling, 
the gas collapses into any low-mass dark matter halos present. 

We model one such halo by using a fixed background p otential. The potential is spher- 
ically symmetric with a density profile (INavarro et al.l 119971 ) 



Pdri 



(^cPcrit 



r/rs(l + r/rs)2' 



(6) 



where is a scale radius, 5c is a characteristic (dimensionless) density and pcrit = SH^/8ttG 
is the critical density for closure. Note that we use a value for the Hu bble constant of 
H = 72kms~^ Mpc~^ (jSpergel et al.ll2003l ). Following iNavarro et al.l (119971 ) we calculate the 
characteristic density and the scale radius from a given redshift and dark halo mass. We 
truncate the halo at the radius at which the value of pdm given by Equation [6] equals the 
cosmological background density at the beginning of the simulation. 

In all of the runs presented here, the halo mass is 7.8 x 10^ Mq and the initial redshift 
is z = 25. We truncate the halo at a radius rt = 0.49 kpc. The virial temperature of the 
resulting halo is 1900 K and the virial radius is 0.1 kpc. In physical units, the scale radius of 
the halo is 29 pc and the full computational volume is a box of side length 1 kpc. 



We use periodic boundary conditions for the hydrodynamic part of the force calculations 
to keep the gas bound within the computational volume. The self-gravity of the gas and the 



- 8 - 



gravitational force exerted by the dark matter potential are not calculated periodically since 
we assume that other dark matter halos and their gas content are distant enough to neglect 
their gravitational influence. 

We begin our simulations with a uniform distribution of gas with an initial density pg, 
taken to be equal to the cosmological background density, and then allow the gas to relax 
isothermally until it reaches hydrostatic equilibrium. This initial phase of the simulation 
is merely a convenient way to generate the appropriate initial conditions for the simulation 
proper, and so we do not include the effects of chemical evolution or cooling during this 
phase. 

The mass of gas present in our simulation was taken to be a fraction f2b/f2dm of the 
total mass of dark matter, where the dark matter density = — Q^,, and where fib 
is the baryon de nsity and i s the matter density. We take values for the cosmological 



parameters from ISpergel et al.l (120031 ) of f2b = 0.047 and f2m = 0.29, giving us a total gas 
mass of Mg = 0.19 Mdm- Mim is the sum of the halo mass and the mass of the dark matter 
background in the simulated volume, and has a value Mdm = 1-84 x 10^ Mq. Therefore, 
Mg = 3.5 X 10^ Mq. Initially, we use 1.4 x 10^ SPH particles to represent the gas, and so 
each particle has a mass rrip = 2.5 Mq. Due to the reflnement, the number of particles rises 
during the simulation up to 9 x 10®, and the SPH particles that represent the gas in the 
collapsing region have a mass of rrip = 0.015 M© (see Appendix |A|) . Our SPH smoothing 
kernel encompasses approximately 50 particles a nd since we need twice this number in order 



to properly resolve gravitationally bound objects (IBate fc Bur kertll 19971 ). our mass resolution 



is Mres ^ 100 mp = 1.5 Mq. 



In the runs presented here we introduce sink particles with an accretion radius of 0.1 pc 
when the density rises above 1.25 x lO^cm"^. At this density Mj exceeds Mj-es as long as 
T > 30 K. As this is much smaller than the CMB temperature at z = 25, we always satisfy 



the iBate &: BurkertI (119971 ) resolution criterion by a comfortable margin. 



We follow the simulations until collapse occurs and a sink forms or as close as possible 
to a Hubble time within reasonable computing time. The Hubble time at a redshift of z = 25 
is ~ 100 Myr. After more than a Hubble time our assumption of an isolated dark matter 
halo will no longer be valid, since most halos will have undergone at least one major merger 
or interaction by this time. The dynamical time within the halo varies between 4 Myr in the 
center and 450 Myr in the outer regions. 
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2.3.2. Metallicity 



We study both zero metallicity gas and gas that has been enriched to 10~^ Zq. A 
metallicity oi Z = 10~^ Zq is an upper limit deri ved from QS Q absorption-line studies of 
the low column density Lyman-a forest at z ~ 3 (|Pettinilll999l ). Estimates of the globally- 
averaged metallicity produced by the sources responsible f or reionization are also typically 
of the order of 10~^ Zq (see e.g. iRicotti fc Ostrikerl |2004| ). We also run a model with a 
metallicity of Z = 10^^ Zq to investigate the influence of an even higher metallicity. 

In our simulations of metal-enriched gas, we assume that mixing is efficient and that 
the metals are spread out uniformly throughout the computational domain. We also assume 
that the relative abundances of the various metals in the enriched gas are the same as in 
solar metallicity gas; given the wide scatter in observational determinations and theoretical 
predictions of abundance ratios in very low metallicity gas, this seems to us to be the most 
conservative assumption. However, variations in the relative abundances of an order of 
magnitude or less will not significantly alter our results. We have denoted runs with zero 
metallicity with "ZO" and runs with metallicity with "Z-3" and "Z-1". 



2.3.3. Rotation 

We further investigate the influence of rotation on the collapse and possible fragmenta- 
tion of the gas. Theoretically, we do expect rotation to become important on scales r -C Xr^ir, 
wh ere ryjr is the virial radius of the halo and A is the dimensionless spin parameter, given 



by (|Peebleslll97lh 



where J is the total angular momentum of the halo, M is the halo mass and E is the total 
(kinetic plus potential) energy of the halo. Previous work has shown that for halos with the 
range of masses and redshifts conside red in this paper, A has a lognormal distribution, with 



mean A = 0.035 (lYoshida et al.ll2003l ). 



In Paper II we showed that rotation has very little effect on the evolution of the gas at 
early times but it did appear to affect the evolution once the density exceeds n = 100 cm~'^, 
significantly slowing the collapse. Here we investigate if the rotation results in the formation 
of a rotationally supported disk. 

We perform three simulations in which the initial gas distribution was given a non-zero 
angular momentum. Within the virial radius of the halo, the gas was placed into rotation 
with constant angular velocity. At larger radii, the initial angular velocity decreased linearly 
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with radius, reaching zero at the truncation radius of the halo. These runs, which we 
designate hereafter using "ROT", had a spin parameter A = 0.05. 



2.3.4- Turbulence 



Finally, we study the influence of low and intermediate levels of turbulence on the evolu- 
tion of the gas within the dark matter halo. Turbulence establishes a network of interacting 
shocks, where converging flows and shear generate filaments of high density. We investi- 
gate if the turbulence can provide high density gas that acts as a seed for fragmentation. 
For one run we assume a turbulent energy of 5% of the internal energy and denote it with 
"TURBl", while for two further runs with turbulence we take a value of 10% and denote 
it with "TURB2". We have included turbulence in our version of the code that is driven 



uniformly with the method described by iMac Low et al.l (119981 ) and iMac Lowl (119991 ). We 
insert energy on scales of the order of the size of our computational domain, i.e. with wave 
numbers k = 1..2. 



2.3.5. Choice of Parameters 

In the choice of parameters, we first restrict ourselves to the primordial case. This 
enables us to directly compare our results with those of Paper IV and to study the influence 
of differences in the initial conditions on the dynamical evolution. Since we suspect that the 
radial symmetry of our initial conditions in run ZO may artificially suppress fragmentation, 
we investigate the effects of adding turbulent energy to the gas: a low level of energy in run 
ZO-TURBl and a higher level in run Z0-TURB2. 

We next investigate the effects of enriching the gas with a low level of metallicity, Z = 
10^^ Zq. In this case, we start our investigation with a run which incorporates turbulence at 
the TURB2 level; since we find no fragmentation in this case, we can reasonably assume that 
we will not find it in runs with less turbulence, or with the spherical symmetry unbroken. 
We therefore do not include the corresponding runs Z-3-TURB1, and Z-3 in our analysis. 

Finally, we study the effects of rotation. Here we consider three different metallicities, 
the purely primordial case (ZO-ROT), as well a.s Z = 10^^ Zq (Z-3=R0T) and Z = 10^^ Zq 
(Z-l-ROT). We add one run with a high metallicity to investigate the influence of the metals 
in the rotating case more thoroughly. We do not mix turbulence and rotation in any of the 
models to prevent confusion of the respective effects. A complete list of runs is given in 
Table [31 
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3. Results 

In Figures [1] and [2] we show the state of the gas after 52 Myr. During this time the gas 
has cooled down from 10^ K, initially mainly due to Lj-a cooling and later due to molecular 
cooling and atomic fine structure cooling where metals are present. In agreement with the 
results presented in Paper IV, the inclusion of molecular hydrogen cooling largely erases 
the difference in temperature structure between zero and low-metallicity gas (see Figured]). 
Only gas with metallicity as high as 0.1 Zq cools at substantially lower density than zero 
metallicity gas, as seen in the rotating models shown in Figure [2l In all cases, the gas is able 
to cool down to the CMB temperature at our initial redshift z = 25. 

In Paper II we found that the centrally concentrated halo profile that we have chosen to 
study suppressed fragmentation at numbe r densities below 500 cm~^ even for metallicities 



far above the threshold Zth suggested by iBromm et al.l (120011 ). Our current models show 
no fragmentation during collapse to number densities as high as 10^ cm~'^ for metallicities 
reaching as high as 10~^ Zq in one rotating case. 

We consider several different variations on the initial conditions of the gas in order to 
understand the lack of fragmentation that we find. Model ZO is a zero metallicity model 
that collapses spherically, with temperatures dropping to the CMB temperature (see Fig. [T]). 
In Paper IV, we showed that fragmentation proceeds readily at these temperatures and 
densities in a top-hat potential with small perturbations, even at zero metallicity, but in the 
centrally concentrated potential studied here, no fragmentation occurs (top panel of Fig. [3]). 
One central sink particle forms and accretes. 



3.1. Turbulence 

We then consider whether breaking spherical symmetry can lead to fragmentation. To do 
this, we initialize the gas with two different levels of turbulent kinetic energy (Table [3]) in runs 
ZO-TURBl and Z0-TURB2. Converging flows effectively produce randomly oriented density 
perturbations, as shown in the middle panels of Figure [31 However, no further gravitational 
fragmentation occurs. In the less turbulent TURBl run, a sink particle forms from collapsing 
gas during the run, while even that does not occur in the more energetic TURB2 runs. 
Increasing the metallicity to Z = 10^^ Zq in run Z-3-TURB2 neither markedly increases the 
cooling compared to the zero metallicity case (Fig. [1]), nor produces fragmentation (bottom 
panel of Fig. [3]) . 
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3.2. Rotation 

Up to this point organized rotation has not been introduced into our models. Disk 
formation naturally leads to fra gmentation and the form ation of multiple stellar systems in 



present-day star formation (e.g. iMatsumoto et al.lll997l ). so we investigate whether it could 



play a similar role in this situation. In Figure H] we show horizontal and vertical number 
density cross sections of runs in this parameter study. We find that thin disk formation 
proceeds only in the Z-l-ROT case, with metallicity Z = 0.1 Zq. In models ZO-ROT and 
Z-3-R0T, with metallicity on either side of the threshold Zth, however, the gas cannot 
cool strongly enough to form a thin disk, and merely forms a slightly flattened, pressure- 
supported structure. In none of these three models does fragmentation set in, nor does the 
central density even grow fast enough for a sink particle to form within 60 Myr, a significant 
fraction of a Hubble time at that redshift. 



4. Discussion 



The basic thermodynamic behavior shown in Figures [T] and [2] underlies our conclusion 
both in Paper IV and here that the introduction of metal line cooling does not yield a 
substantial difference in star formation behavior once H2 and HP cooling is als o taken into 
account. Recall that the fragmentation threshold found by lBromm et al.l (1200 ll ) occurred in 
simulations neglecting these coolants. 

In the centrally condensed potential that we study here, a wide variety of initial condi- 
tions for the gas yield a monolithic central collapse. Figures |5] and [6] show that the central 
density depends on the initial conditions. Turbulence slows the growth of the central density 
slightly (Fig. [5]), but both spherically symmetric and turbulent initial conditions collapse and 
form a sink particle. The mass in dense gas (over an arbitrary value of nthres = 10^ cm"^) 
grows rapidly in both cases once collapse sets in (Fig. \5^. 

Rotation allows for the formation of stable disks over large fractions of the local Hubble 
time. Figure [6] shows that once a rotationally supported structure has formed, the central 
density and amount of mass at high density ceases to grow quickly. This is because the 
continued growth of the central density now depends on angular momentum transport within 
the rotating structure. Molecular viscosity is sufficiently low that no further accretion would 
happen at all if that were the only cause of transport. In modern star formation, the 
additional transport in circumstellar disks occurs be cause of magne torotational instability 
(IBalbus fc Hawleylll998l ) or gravitational instability (jGammiell200ll ). However, our models 
include no magnetic fields, so the magnetorotational instability cannot act here. 
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4.1. Gravitational Instability 

To examine whether gravitational inst ability rnay be important in driving accretion or 



fragmentation, we examine the value of the iToomrd (Il964j ) instability parameter 



Q = kCs/ttGS, (8) 

where G is Newton's gravitation constant, S is the disk surface density, k is the epicyclic 
frequency, and Cg is the sound speed. Radial gravitational instability occurs for Q < 1. In 
Figure [7] we examine our rotating runs. We see that, in agreement with the observed lack 
of fragmentation, the disks are strongly Toomre stable for the first 50 Myr, with Q > 10 
throughout, and only begin to reach the unstable regime at the very center thereafter. 

In fact, given the lack of gravitational instability or magnetic fields, the accretion towards 
the center seen in the simulation still represents only an upper limit on the actual accretion 
rate, as all the angular momentum transfer occurring must be due to numerical viscosity. In 
reality, virtually all galaxies interact with other galaxies over the course of a Hubble time. 
The resulting tidal perturbations drive strong accretion and subsequent star formation. Our 
results suggest that, at least at early times when only small halos are collapsing, galaxy 
interactions may be essential for driving star formation. 

Another way of looking at the development of gravitational instability is to examine 
the total number of Jeans masses at every number density. We space our data equally in 
logarithmic number density with a bin size of 0.14 and in each bin we calculate the Jeans 
mass 

M, = |po(|)' = ^G-/V„-"V, (9) 

where we define the Jeans mass Mj as the mass originally contained within a sphere of diam- 
eter A J and uniform density po- -^J denotes the critical wavelength for the Jeans instability. 
Since the sound speed Cg depends on the temperature and chemical composition, we use the 
mean value of Cg in each density bin. We calculate the number of Jeans masses A^j above a 
certain density uq as 

Nj = M{n>no)/M,{no). (10) 

We show the result in Figure [SI This compares the standard primordial run ZO to the highest 
metallicity rotating run Z-l-ROT. Whereas ZO reaches Jeans number unity at the highest 
densities at t > 55 Myr, Z-l-ROT only approaches Jeans number unity at lower densities 
n ~ 100 cm~'^. This gas is widely distributed away from the center of the disk, however, so 
it cannot collect and collapse towards a common center. Too little material reaches the high 
densities at the very center of the disk to start runaway collapse there. 
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4.2. Metallicity Thresholds 



How does our result agree with oth er work on low-met a llicity star formation? We have 
already discussed the comparison with iBromm et al.l (120011 ) . The same authors examined 
their zero metallicity case in the absence of an ambient UV radiation sufficient to dissociate 
H2, so that molecules can form and contribute to the cooling. They found the same result as 
was later reached b y Paper IV, that th e top hat potential does lead to fragmentation even 
at zero metallicity (IBromm et al.ll2002l ). 



A more advanced computation was performed bvlSmith fc Sigurdssonl (120071 ). who used 



the adaptive mesh grid code Enzo (lO'Shea et al.l 12004 ) and a comprehensive treatment of 



metal cooling rates as detailed by ISmith et al.l (j2008aj). Instead of a top-hat dark matter 
potential, they used cosmological i nitial conditions, with collapse followed over 28 levels 
of refinement. Like iBromm et al.l ((2002) they assumed no dissociating radiation and full 
molecular cooling. As a result of their choice of initial conditions, their collapsing objects 
were centrally concentrated, similar to ours. In agreement with our results, they found no 
fragmentation below densities of 10^ cm~'^. At higher densities of 10^ cm~^, they did see hints 
of fragmentation of their central object into two to four objects at metallicities Z > 10~'^ Zq, 
but that alone seems insufficient to mark a shift from a primordial to a modern initial mass 
function. 



Following up on that work ISmith et al.l (l2008bl ) presented a suite of simulations inves- 
tigating the critical metallicity threshold. Two of their runs show fragmentation below den- 
sities of 10^ cm~^ into 2 bound clumps but most of their runs show fragmentation at higher 
densities. Their set of models which collapses at a redshift similar to the one used in our 
calculations only results in 2 fragments at a metallicity of Z = 10~^'^ Zq and Z = 10~^ Zq. 
These differences could also be due to the different initial conditions used. Our simulations 
began with hot, ionized gas, whereas their simulations started with cold neutral gas. In 
Pape r II we explained th e consequences of these differences in more detail. In agreement 
with ISmith et al.l (l2008bl ) we also find that the gas reaching the CMB temperature results 
in a thermally stable gas cloud. 



Schneider et al.l (120021 . |2006| ) and lOmukai et al.l (120051 ) propose that rather than atomic 



or molecular line cooling determining the metallicity where strong fragmentation sets in, 
dust continuum cooling does so. Surprisingly, their estimates of cooling rates suggest that 
fragmentation at high densi t ies n > 10^° cm~^ already sets in for metallicities as low as 



Z > 10 ^ Zq. lOmukai et al.l (120051 ) demonstrates that even at these low metallicities, dust 



cooling gives an effective adiabatic index under unity for a range of densities 10 < n < 
10^^ cm~^, promoting rapid fragmentation of objects collapsing through that density range. 
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High resolution nested-grid models by lMachida et al.l (12008) show bi n aries forming even 



at primordial metallicity Z = ai densities n > 10^** cm""^. iMachidal (120081 ) then found 
that binary fragmentation is actually enhanced at decreasing metallicity, with the rotation 
threshold for fragmentation aX Z = lying a full two orders of magnitude lower than for 
solar metallicity Z = Zq. Low metallicity fragmentation occurs at densities as high as 
10^^ cm~^ for Z = 0, in contrast to modern star formation where binaries only form for 
10^^ < n < 10^^ cm~^. However, primordial binary formation still yields very massive stars, 
simply in pairs rather than alone. 



Further simulation of collapsing cores using the lOmukai et al.l (120051 ) equation of state 
has shown that fragmentation promoted by dust cooling at high densities indeed appears 
able t o shift star formatio i i frorn a primordial to a more present-day-like initial mass func- 
Tsuribe &: Omukail (120061 ) followed fragmentation of low metallicity Z < 10^'^ Zq, 



tion. 



high density n > 10 cm gas in the non-rotating case. IClark et al.l (120081 ) included rota- 



tion, and used sink p articles to follow the collapse of a complete cluster. They agree with 
Machida et al.l ( 20081 ) in finding only small numbers of high-mass fragments in the primordial 
case, but already at metallicities of Z ~ 10~^ Zq, they find an average mass under 1 Mq 
approaching a modern initial mass function. 



0! 



Dust cooling looks likely to be the primary physics determining the shift in typical stellar 
mass from primordial to modern stars, with a threshold of Z ~ 10~^ Zq consistent with the 
lowest metallicity stars observed in the modern galaxy. Conversely, as we have shown in 
Paper IV and this paper, atomic line cooling does not appear to dominate this transition, 
and it appears unlikely that the transition occurs at metallicities as high a.s Z = 10~^ Zq. 
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A. Particle splitting 

In Paper II we presented low-resolution simulations that showed that the density in- 
creases most rapidly close to the center of the dark matter halo. On-the-fiy splitting 



(IKitsionas fc WhitworthI |2002| ) with two levels of refinement at 2 different radii provides 
us with the highest mass resolution at the region of interest. Every SPH particle that 
crosses a radius ri of 0.3 kpc towards smaller radii acts as a parent particle and spawns 13 
child particles with a mass of mp/13. The same happens at a radius r2 of 0.2 kpc where 
gas particles crossing r2 are split again. The mass of the gas particles within r2 is thus 
mp/169 ~ 1.5 X 1O~^M0. In the region of collapse we can thus resolve Jeans masses down 
to 1.5 Mq, as compared to Mj-es = 200 Mq in the low resolution simulations. Each child 
particle inherits directly velocity, internal energy and fractional abundances of the chemical 
species. Subsequently the child particles are evolved with standard SPH procedures. The 
only difference is that, to mitigate interactions between adjacent particles having different 
masses, we have modified the scheme by which GADGET calculates the smoothing lengths 
of particles. We now evolve hi, the smoothing length for particle i, so that the radius hi 
contains between Ancigh — Adev and Aneigh + Adev, where Adev is the allowed deviation in 
the number of neighbors, as well as a mass between Ancigh — Adev and Aneigh + Adev times 
the mass of particle i. If both conditions cannot be satisfied simultaneously we allow for a 
larger number of neighbors or a larger amount of mass within a sphere of radius hi. For one 
specific set of parameters we carried out two simulations, one with the refinement technique 
and the other without it, and verified that temporal perturbations caused by refinement do 
not affect the results presented in Section [31 
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Fig. 1. — Gas temperature vs. number density for runs ZO {red solid line), ZO-TURBl [blue 
dotted line), Z0-TURB2 [green dashed line), and Z-3-TURB2 [black dot-dashed line). The 
thin hues show the la-deviation around the mean. The time of all plots is t = 52 Myr, 
except for plot Z0-TURB2, which is at t = 50 Myr. At this time the gas has cooled down 
from the initial temperature of 10^ K, initially mainly due to Lj-a cooling and later also 
due to molecular cooling and atomic fine structure cooling if metals are present. We show 
the CMB temperature with the horizontal dashed line. Runs ZO and ZO-TURBl reach the 
CMB temperature in t = 52 Myr just from H2 cooling despite their lack of metals. On the 
other hand, run Z-3-TURB2, which has a higher fraction of turbulent energy, does not reach 
the CMB temperature in t = 52 Myr despite having metallicity Z= 10~^ Zq. 
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Fig. 2. — Gas temperature vs. number density for runs ZO-ROT {red solid line), Z-3-R0T 
{blue dotted line), and Z-l-ROT {green dashed line). The thin hues show the lex- deviation. 
The halos with rotation have a spin parameter of 0.05. The time of the plots is t = 52 Myr 
as in Figure [H We show the CMB temperature with the horizontal dashed line. The run 
Z-l-ROT reaches the CMB temperature at a relatively low density value due to the high 
content of metals that contribute effectively to the cooling. The other two runs show only 
little difference at higher densities due to the metallicity of run Z-3-R0T. 
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Fig. 3. — From top to bottom: Cross sections showing number density for runs ZO, ZO- 
TURBl, Z0-TURB2, and Z-3-TURB2 at 50 Myr. Left: x-y plane, cut at z = 0.5 kpc. Right: 
x-z plane, cut at y = 0.5 kpc. 



0.49 0.495 0.5 0.505 0.49 0.495 0.5 0.505 
X [kpc] X [kpc] 




2 4 2 4 

og density [cm~'^] log density [cnn""^] 



Fig. 4. — From top to bottom: Cross sections showing number density for runs ZO-ROT, 
Z-3-R0T, and Z-l-ROT, at 50 Myr. Left: x-y plane, cut at ^ = 0.5 kpc. Face-on view. 
Right: x-z plane, cut at y = 0.5 kpc. Edge-on view. 
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Fig. 5. — Left panel: Maximum gas number density vs time for runs ZO {solid line), ZO- 
TURBl {dotted line), Z0-TURB2 {dashed line), and Z-3-TURB2 {dot-dashed line). Note 
that the hnes for runs Z0-TURB2 and Z-3-TURB2 overlap completely up to t = 50 Myr. 
Right panel: Mass over density threshold nthres = lO^cm"^. 




Fig. 6. — Left panel: Maximum gas number density vs time for runs ZO-ROT {solid line), 
Z-3-R0T {dotted line), and Z-l-ROT {dashed line). Right panel: Mass over density threshold 
nthres = lO^cm"^. 



- 25 - 



r 



5 Myr 
15 Myr 
25 Myr 
35 Myr 
45 Myr 
55 Myr 
65 Myr 
86 Myr 



f V II ill 



I i 



^ — I — I — I — I — I — I — I — I — — I — I — I — I — I — I — I- 




-I — I — I — I — I — I — I — I— 



"■■is 



■J — I — I — I — I- 




I s 



>r: I . 



5 Myr 
15 Myr 
25 Myr 
35 Myr 
45 Myr 
55 Myr 
65 Myr 
82 Myr 



li 




-4 -3 -2 

logarithmic radius [kpc] 



Fig. 7. — From top to bottom: Radial distribution of the Toomre parameter Q for the runs 
with rotation ZO-ROT, Z-3-R0T, and Z-l-ROT. The value of Q is shown on a logj^g scale. 
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Fig. 8. — We show the temporal evolution of the number of Jeans masses as a function of 
number density for the runs ZO {solid lines) and Z-l-ROT {dashed lines). The runs Z-l-ROT 
reach a Jeans number of 1 {dot- dashed horizontal line) but not at the high number densities 
that are required for efficient cooling and collapse. In contrast the runs ZO show a clear 
trend of the Jeans number rising with time at high densities. The data are equally spaced 
in logarithmic number density with a bin size of 0.14. 
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Table 1. 



Species Method of solution 



e 


Conservation law 


H+ 


Rate equation 


H 


Conservation law 


H- 


Equilibrium abundance 


H+ 


Equilibrium abundance 


H2 


Rate equation 




Equilibrium abundance 


D+ 


Rate equation 


D 


Conservation law 


HD 


Rate equation 


He+ 


Rate equation 


He 


Conservation law 


C+ 


Rate equation 


c 


Conservation law 


c- 


Equilibrium abundance 


0+ 


Rate equation 





Conservation law 


0- 


Equilibrium abundance 


Si+ 


Rate equation 


Si 


Conservation law 


Si++ 


Rate equation 


CH+ 


Equilibrium abundance 


CH+ 


Equilibrium abundance 


CH3+ 


Rate equation 


CH 


Rate equation 


CH2 


Rate equation 


0H+ 


Equilibrium abundance 


H2O+ 


Equilibrium abundance 


H3O+ 


Equilibrium abundance 


OH 


Rate equation 


H2O 


Rate equation 
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Table 1 — Continued 



Species Method of solution 



C0+ 


Equilibrium abundance 


CO 


Rate equation 


HOC+ 


Equilibrium abundance 


HCO+ 


Rate equation 


0+ 


Equilibrium abundance 


O2 


Rate equation 


C2 


Rate equation 


SiH+ 


Equilibrium abundance 



Table 2. Processes included in our thermal model. 



Process 



References 



Cooling: 

Lyman-Q cooling 

He electronic excitation 

Thermal bremsstrahlung 

Compton cooling 

H2 rotational, vibrational lines 

HD rotational, vibrational lines 

Fine structure lines (C, C"*", O, Si, Si^) 

CO rotational, vibrational lines 

H2O rotational, vibrational lines 

OH rotational, vibrational lines lines 

H"*" recombination 

He"*" recombination 

H &; He collisional ionization 

H2 collisional dissociation 

Heating: 

H2 gas-phase formation 



Cen (1992) 
Cen (1992); Brav et al. (2000 ^ 
Shapiro k Kang (1987 '! 
Cen (1992) 

Le Bourlot et al. (1999); Ripamonti k Abel (2004) 
Lipovka. Nuhez-Lopez. &: Avila-Reese (2005) 
Many sources; see papers I, III 
Neufeld k Kaufman (1993 ): Neufeld. Lepp k Melnick (1995 ) 
Neufeld k Kaufman (1993 ): Neufeld. Lepp k Melnick (1995) 
Pavlovski et al. (2002) 
Ferland et al. (1992) 
Aldrovandi k Peauignot (1973 ) ; Hummer k Storev (1998) 
Janev et al. (1987) 
Many sources; see papers I, III 

Many sources; see papers I, III 
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Table 3. Physical state of the densest gas within the scale radius rg at time tend 



Run 






-E'turb'^ 


4- d 


rp C 

c,min 


f 

^c,max 




2,C, 


,max 


g 




(Ze) 






(Myrs) 


(K) 


(cm-3) 










ZO 


0.0 


0.0 


0.0 


53 


79 


^sink 


2.9 


X 


10- 


-3 


ZO-TURBl 


0.0 


0.0 


0.05 


71 


66 


^sink 


3.2 


X 


10- 


-3 


Z0-TURB2 


0.0 


0.0 


0.1 


50 


82 


3.5 X 10^ 


2.8 


X 


10- 


-3 


Z-3-TURB2 


10-3 


0.0 


0.1 


65 


66 


1.6 X 10^ 


3.0 


X 


10- 


-3 


ZO-ROT 


0.0 


0.05 


0.0 


86 


67 


2.5 X 10^ 


2.9 


X 


10- 


-3 


Z-3-R0T 


10-3 


0.05 


0.0 


70 


59 


1.5 X 10^ 


2.9 


X 


10- 


-3 


Z-l-ROT 


10-^ 


0.05 


0.0 


82 


67 


1.9 X 10"^ 


2.7 


X 


10" 


-3 



'^Metallicity of the gas. 
''Spin parameter of the halo. 

Added turbulent energy in units of the initial internal energy. 
•^Time at the end of the simulation. 

^Minimum temperature of the gas within the scale radius Tg. 

^Maximum number density of the gas within the scale radius rg. rigink = 
1.25 X 10^ cm-3 jg ^jjg threshold density for sink creation. 

^Maximum fractional H2 abundance within the scale radius Tg. 



